A CONVERGENT LINEAR FINITE ELEMENT SCHEME FOR THE 
MAXWELL-LAND AU-LIFSHITZ-GILBERT EQUATION 



L\ BANAS, M. PAGE, AND D. PRAETORIUS 

Abstract. We consider a lowest-order finite element discretization of the nonlinear system 
of Maxwell's and Landau-Lifshitz-Gilbert equations (MLLG). Two algorithms are proposed 
to numerically solve this problem, both of which only require the solution of at most two 
linear systems per timestep. One of the algorithms is fully decoupled in the sense that each 
timcstcp consists of the sequential computation of the magnetization and afterwards the 
magnetic and electric field. Under some mild assumptions on the effective field, we show 
that both algorithms converge towards weak solutions of the MLLG system. Numerical 
experiments for a micromagnetic benchmark problem demonstrate the performance of the 
proposed algorithms. 



I. Introduction 

The understanding of magnetization dynamics, especially on a microscale, is of utter rele- 
vance, for example in the development of magnetic sensors, recording heads, and magneto- 
resistive storage devices. In the literature, a well accepted model for micromagnetic phe- 
nomena, is the Landau-Lifshitz-Gilbert equation (LLG), see (la). This nonlinear partial 
differential equation describes the behaviour of the magnetization of some ferromagnetic 
body under the influence of a so-called effective field. Existence (and non-uniqueness) of 
weak solutions of LLG goes back to [AS, V'85]. Existence of weak solutions for MLLG was 
first shown in [CF'98]. For a complete review of the analysis for LLG, we refer to [C, GC, KP] 
or the monographs [HS, P] and the references therein. As far as numerical simulation is con- 
cerned, convergent integrators can be found e.g. in the works [BP, BKP] or [BBP'08], where 
the latter considers a weak integrator for the coupled MLLG system. From the viewpoint 
of numerical analysis, the integrator from [BKP] suffers from explicit time stepping, since 
this imposes a strong coupling of the timestep-size k and the spatial mesh-size h. The in- 
tegrators of [BP, BBP'08], on the other hand, rely on the implicit midpoint rule for time 
discretization, and unconditional convergence is proved. In practice though, a nonlinear sys- 
tem of equations has to be solved in each timestep, and to that end, a fixed-point iteration 
is proposed in the works [BP, BBP'08]. This, however, again leads to a coupling of h and k, 
and thus destroys unconditional convergence. Using the Newton method for the midpoint 
scheme seems to allow for larger time-steps, see [DSM] and also [BBrP'12], where an efficient 
Newton-multigrid nonlinear solver has been proposed. 

In [A'08], an unconditionally convergent projection-type integrator is proposed, which, 
despite the nonlinearity of LLG, only requires the solution of one linear system per timestep. 
The effective field in this work, however, only covers microcrystalline exchange effects and 
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is thus quite restricted. In the subsequent works [A'll, GHMPS, GPS] the analysis for this 
integrator was widened to cover more general (linear) field contributions, where only the 
highest-order exchange contribution is treated implicitly, whereas the other contributions 
are treated explicitly. This allows to minimize computational effort while still maintaing 
unconditional convergence. Finally, in the very recent work [BSFFGPP'12], the authors 
could show unconditional convergence of this integrator, where the effective field consists of 
some general energy contributions, which are only supposed to fulfill a certain set of prop- 
erties. This particularly covers some nonlinear contributions, as well as certain multiscale 
problems. In addition, it is shown in [BSFFGPP'12] that errors arising due to approximate 
computation of field contributions like e.g. the demagnetizing field can be incorporated into 
the analysis. 

In [A'12], the authors also investigate a higher-order extension of this algorithm which, 
however, requires implicit treatment of nonlocal contributions like the magnetostatic stray- 
field. 

In our work, we extend the analysis of the aforementioned works and show that the integra- 
tor from [A'08] can be coupled with a weak formulation of the full Maxwell system (lb)-(lc). 
For the integration of this system, we propose two algorithms that only require the solution 
of one (Algorithm 2) resp. two linear systems (Algorithm 3) per timestep while still guar- 
anteeing unconditional convergence (Theorem 6). The contribution of the present work can 
be summarized as follows: 

• We extend the linear integrator from [A'08, GHMPS] to time-dependent contribu- 
tions of the effective field by considering the full Maxwell equations instead of the 
magnetostatic simplification. 

• Unlike [BBP'08], at most two linear systems per timestep, instead of a coupled nonlin- 
ear system, need to be solved. Nevertheless, we still prove unconditional convergence. 

• Unlike [BBP'08], the decoupling of the Maxwell and the LLG part in the integrator of 
Algorithm 3 is rigorously included into the convergence analysis of the time-marching 
scheme. 

Outline. The remainder of this paper is organized as follows: In Section 2, we recall 
the mathematical model for the full Maxwell-LLG system (MLLG) and recall the notion of 
a weak solution (Definition 1). In Section 3, we collect some notation and preliminaries, 
as well as the definition of the discrete ansatz spaces and their corresponding interpolation 
operators. In Section 4, we propose two algorithms (Algorithm 2 and 3) to approximate the 
MLLG system numerically. The large Section 5 is then devoted to our main convergence 
result (Theorem 6) and its proof. Finally, in Section 6, some numerical results conclude this 
work. 



2. Model Problem 

We consider the Maxwell-Landau-Lifshitz-Gilbert equation (MLLG) which describes the 
evolution of the magnetization of a ferromagnetic body that occupies the domain u d Q C 
IR 3 . For a given damping parameter a > 0, the magnetization m:(0,T)xu-f§ 2 and the 
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electric and magnetic fields E, H : (0, T) x ft ->■ IR 3 satisfy the MLLG system 



m t — am x m ( = — m x H e g in ojt '■= (0, T)xw (la) 
e E t - V x H + crx„E = - J in ft T := (0, T) x ft (lb) 
/i H t + V x E = — /ionij in ft T , (lc) 

where the effective field H c g- consists of H c g = C e Am + H + 7r(m) for some general energy 
contribution 7r which is assumed to fulfill a certain set of properties, see (13)-(14). This 
is in analogy to [BSFFGPP'12]. We stress that, with the techniques from [BSFFGPPT2], 
an approximation ix^ of n can be included into the analysis, as well. We emphasize that 
throughout this work, the case H cff = C e Am + H + C a D$(m) + H ext is particularly covered. 
Here, <£>(•) denotes the crystalline anisotropy density and H ex t is a given applied field. The 
constants Eq, /io > denote the electric and magnetic permeability of free space, respectively, 
and the constant a > stands for the conductivity of the ferromagnetic domain u. The field 
J : ft r — > IR 3 describes an applied current density and Xu> '■ ^ — > {0, 1} is the characteristic 
function of u. As is usually done for simplicity, we assume ft C IR 3 to be bounded with 
perfectly conducting outer surface <9ft into which the ferromagnet u (s ft is embedded, and 
fl\cu is assumed to be vacuum. In addition, the MLLG system (1) is supplemented by initial 
conditions 

m(0, •) =m° inw and E(0, •) = E°, H(0, •) = H° in ft (Id) 

as well as boundary conditions 

<9 n m = on du T , E x n = on <9ft T . (le) 

Note that the side constraint |m| = 1 a.e. in ojt does not need to be enforced explicitely, 
but follows from |m°| = 1 a.e. in ui and <9 t |m| 2 = 2m • m t = in oj t , which is a consequence 
of (la). This behaviour should also be reflected by the numerical integrator. In analogy to 
[CF'98, BBP'08], we assume the given data to satisfy 

m° G H l (u, § 2 ), H°, E° G L 2 (ft, IR 3 ), JGL 2 (ft T ,R 3 ) (If) 

as well as 

div(H° + Xcum ) =0 in ft, (H° + x,m°, n) = on <9ft. (lg) 
With the space 

H (curl,ft) := G L 2 (ft) : V x up G L 2 (ft),v? x n = on T}, 
we now recall the notion of a weak solution of (la)-(lc) from [CF'98]. 

Definition 1. Given (If) — (lg), the tupel (m, E, H) is called a weak solution of MLLG if, 
(i) m G H^wy) with |m| = 1 almost everywhere in ojt and (E, H) G L 2 (Qt); 
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(ii) for all if e C°°{oo T ) and ( e C c °°([0, T); C°°(Q) n H°(curl, SI)), we have 

/ (m t ,(p)-a ((m x m t ),<p) = -C e / ((Vmxm),V^) (2) 



+ ((H x m).ip) + / ((i(m)xm),^ 
-e I (E,Ct>-/ (H,Vx() + J (E,C> = -/ (J,C>+eo/<E°,C(O,0>, (3) 

-/'«» /' (H,C*>+/ (E.VxC) = -//u /' (m t ,C)+^o /(H°,C(0,.)>; (4) 

(iii) t/iere /io/ds m(0, •) = m° m £/ie sense o/ traces; 

(iv) /or almost all t' G (0,T), we /iawe bounded energy 

||Vm(f + ||m t ||2 2K) + ||H(0||L. (n) + ||E(OllL. ( n) < C (5) 

where C > is independent oft. 

Existence of weak solutions was first shown in [CF'98]. We note, however, that our analysis 
is constructive in the sense that it also proves existence. 

Remark. Under additional assumptions on the general contribution 7r(-) ; namely that ir(-) 
is self-adjoint with \\n(n) ||l 4 (w) < C / or n £ L 2 (w) itrai/i |n| < 1 almost everywhere, the 
energy estimate (5) can be improved. The same techniques as in [BSFFGPP'12, Lemma A.l] 
then show for almost all t' G (0, T) and e > 

£(m,H,E)(0 + 2(a-eVo||n^ (J,E), 

tt>/jere 

£(m, H, E) := fi C e || Vm|| 2 2(w) + /x ||H||£ 2(n) + £o||E||£ 2(aj) - /i (vr(m), m). 

This is in analogy to [BBP'08]. In particular, the above assumptions are fulfilled in case of 
vanishing applied field H ex t = and ifir(-) denotes the uniaxial anisotropy density. □ 



3. Preliminaries 

For time discretization, we impose a uniform partition = to<^i<---<^v = ^of the 
time interval [0, T\. The timestep size is denoted by k = kj := tj+i — tj for j — 0, . . . , N — 1. 
For each (discrete) function <p, tpi := <p(tj) denotes the evaluation at time tj. Furthermore, 
we write d,(p-' +] := (<p j+1 — <p j )/k for j > 1, and ipi +1 l 2 := (<p j+1 +(p j )/2 for j > and a 
sequence {^}j>o- 

For the spatial discretization, let T h n be a regular triangulation of the polyhedral bounded 
Lipschitz domain Q C IR 3 into compact and non-degenerate tetrahedra. By Th, we denote 
its restriction to w H fi, where we assume that u is resolved, i.e. 

T ft = 7^| w = {T G 7^ : Tflw^f)} and w= (J T. 
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By S l (Th) we denote the standard T^-FEM space of globally continuous and piecewise afflne 
functions from u to R 3 

S\%) := {<f> h G C(U,R 3 ) : </> h \ K G Vi(K) for all K G %}. 

By Xh : C(Q) — > 5 1 (7/ l ), we denote the nodal interpolation operator onto this space. Now, let 
the set of nodes of the triangulation 7^ be denoted by Mh- For discretization of the magne- 
tization m in the LLG equation (la), we define the set of admissible discrete magnetizations 
by 

M h := {<j> h G S\T h ) : \4> h {z)\ = 1 for all z G N h }. 

Due to the modulus constraint |m(t)| = 1, and therefore m t ■ m = almost everywhere in 
ut, we discretize the time derivative v(t,) := m t (t,) in the discrete tangent space which is 
defined by 

K+ h := {j> h G S\T h U : Vft(z) • <Mz) = for all z G A4} 

for any <j) h G .M^. For two vectors x, y G 1R 3 , x • y stands for the usual scalar product in M 3 . 

To discretize Maxwell's equations (lb)-(lc), we use conforming ansatz spaces X h C 
H°(curl;fi), y h C L 2 (fi) subordinate to T h n which additionally fullfil Vx^c y h . In 
analogy to [BBP'08], we choose first order edge elements 

X h := {<p h G H°(curl; Q) : tp h \ K G ^i(X) for all X G 7^} 

and piecewise constants 

y h := {(h G L 2 (fi) : ( h \ K G P o (*0 for all K G 7^}, 

cf. [M'03, Chapter 8.5]. Associated with X h) let Z Xfl : H 2 (f2) — )• denote the correspond- 
ing nodal FEM interpolator. Moreover, let 

Xy h : L 2 (fi) y h 

denote the L 2 -orthogonal projection, characterized by 

(C - Iy h C y h ) = for all C G L 2 (Q) and y h G y h . 
By standard estimates, see e.g. [M'03, BS], one derives the approximation properties 

\\<P - Zx h <p\\v>(n) + h\\V x (<p - I Xh <p)\\L,>(n) < C h 2 \\V 2 <p\\ L 2 (n) (6) 
\\C-Zy h C\\v>V)<Ch\\C\\ ma) (7) 

for all tp G H 2 (fT) and ( G H 1 ^). 

Finally, given two expressions A and _B, we write A < £> if there exists a constant c > 
which is independent of h and fc, such that A < cB. In the case A < B and 5 < A, we 
write A ~ B. 

4. Numerical algorithms 

We recall that the LLG equation (la) can equivalently be stated by 

am t + m x m ( = H cff - (m • H cff )m (8) 

under the constraint |m| = 1 almost everywhere in Qt- This formulation will now be used 
to construct the upcoming numerical schemes. Following the approaches of Alouges et 
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AL. [A'08, A'll] and Bruckner et al. from [BSFFGPP'12], we propose two algorithms 
for the numerical integration of MLLG, where the first one follows the lines of [BBP'08]. 

4.1. MLLG integrators. For ease of presentation, we assume that the applied field J is 
continuous in time, i.e. J G C([0, T]; L 2 (fi)) so that «F := J (£_,-) is meaningful. We emphasize, 
however, that this is not necessary for our convergence analysis. 

Algorithm 2. Input: Initital data m°, E° ; and H° 7 parameter 9 G [0,1], counter 
j = 0. For all j = 0, . . . , N — 1 iterate: 
(i) Compute unique solution (v J h , E{ +1 , H{ +1 ) G (/C m j , X h , 34) such that for all {4>h,i>hXh) £ 
x <*ft x yh holds 



«K,W + ((mj x vi) A) = -C7 e (VK + 0kv J h ),V<j> h ) 

+ (Hi +1/2 ,^) + (7r(mi),^), 



(9a) 



coC^Ef 1 ,^) - (Hi +1/2 , V x ^) + ff (x w Ei + 1/2 M = -(J j+1/2 ,^ h ), (9b) 

/i (^Hi +i ,a) + (v x Ei +i/2 ,a) = -//oK,a). (9c) 

m*^ fzl i jtv^ (z^ 

(ii) Define m{ +1 G A4 fe nodewise by m{ +1 (z) = — h -^- for all z E Af h 

|mi(z) + kv{(z)\ 

For the sake of computational and implementational ease, LLG and Maxwell's equations 
can be decoupled which leads to only two linear systems per timestep. This modification is 
explicitely stated in the second algorithm. 

Algorithm 3. Input: Initital data m°, E° 7 and H° ; parameter 9 G [0,1], counter 
j = 0. For all j = 0, . . . , N — 1 iterate: 

(i) Compute unique solution vj t G /C mJ such that for all (j)h G /C mJ holds 

«KA) + (K x v£)A) = -C e (V(m{ + V* fc ) 

(10a) 

+ (ff^ k )+( 1 r(mjl) 1 ^). 

(ii) Compute unique solution (E{ +1 , Hi +1 ) G (A^, 34) snc/j that for all (^hXh) 34 

5oKEf \^) - (Hf 1 , V x V.) + <r( Xw Ei +1 ,^) = -(J'> fc ), (10b) 
M^Hi +1 ,O0 + (V x E£ + \C*) = -^oK,^). (10c) 

nx - ^ fzl | jtv^ (zl 

(iii) De/me m{ +1 G A^ ft nodewise by m{ +1 (z) = — for a// z G A4- 

* K(z) + K(z)| 

4.2. Unique solvability. In this brief section, we show that the two above algorithms are 
indeed well defined and admit unique solutions in each step of the iterative loop. We start 
with Algorithm 2. 

Lemma 4. Algorithm 2 is well defined in the sense that in each step j = 0, . . . , N — 1 of the 
loop, there exist unique solutions (m{ +1 , v{, E^ +1 , H^ +1 ). 
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Proof. We multiply the first equation of (9) by no and the second and third equation by 
some free parameter C\ > to define the bilinear form a J (-, •) on (/C j , X h , 34) by 

h 

a»" ((*,*,©), (MO) ■= aiMt(*,4>)+iM> ((mj x + ^C e 9k (V*, V0) - ^(*,C) 

+ ^i)-^,VxC) + ^i) 

+ ^pce.o + y(v x e.</o + d/io (e,0) 

and the linear functional LP{-) on (/C j , A^, 34) by 

ft 

L J '((0,V,O) := -/ioC e (Vmi,V0) + ^(H£,0) +/i (7r(mJ),*) 

- ft^^) - ^(E^) + ^( H J, V x ^) - ^(xA*) 

-^f (Hi,C)-^(VxEi, . 

To ease the readability, the respective first lines of these definitions stem from (9a), the 
second from (9b), and the third from (9c). Clearly, (9) is equivalent to 

a'"((vi,E£ +1 ,Hj +1 ), (<f> h ,i/, h ,Ch)) = L((<f> h ,i> h ,Ch)) for all {4>hAhXh) e /C m , x X h x y h . 

Next, we aim to show that the bilinear form a J "(-, •) is positive definite on /C j x ^ x ^. 

ft. 

Usage of the Holder inequality reveals that for all (tp,ipX) e x ^ x 34 it holds that 

h 

a* ({MX), MX)) = <*IM>(*,+) + tt>(K x d>),d>)+^C e 9k(Vd>,Vd>) - ^ (0,C) 



+ ^(<M) - ylCVx V) + ^(*M) 
+ ^(CX) + y(V x V,C) + C^o (0,0 
= c^ (M +noC e 9k(y4>, V0) + (d/io - ^)(0,O 

+ ^(^ + %WV0 + ^f(C,O 

> fa-2e(C 1 -l/2))^ ||^||i 2M + ^ll^lli,^ + (f - /io 

=:a , 



where we have used 



(^0>-2£Hli 2H -^IICIIi 2M - 



In order to prove the desired result, we have to show a, b > by choice of C±. We make the 
following ansatz: Fix C2 > k > 0, choose e > in such a way that < aC^ — 2e 2 — 2ae, and 
define C\ := 62/(262 — 4e). Then, Condition a > is equivalent to 

/ 1 Co 1 \ Co — Co -\- 2,£ o 

0<fl = a-2e(-^-^--) =a-e ^ _ 2g <=► < aC 2 - 2ae - 2e 2 



which is true due to the choice of e. The condition b > now gives 

0<6=^-^l k{C l -\)<2eC l ^ k<^\ = C 2 , 

which is automatically fulfilled due to the original choice of k < C 2 - In particular, (9) thus 
admits a unique solution (v{, E^ +1 , H^ +1 ) in each step of the loop. From e K j and the 

Pythagoras theorem, we get |m^(z) + A;v^(z)| 2 = |m^(z)| 2 + A;|v^(z)| 2 > 1. Hence, also step 
{%%) of Algorithm 2 is well defined. This concludes the proof. □ 

The following lemma states an analogous result for the second algorithm. 

Lemma 5. Algorithm 3 is well defined in the sense that it admits a unique solution at each 
step j — 0, . . . , N — 1 of the iterative loop. 

Proof. For the first equation (10a), we define the bilinearform a^(-,-) : K, j x K, j — >■ R by 

h h 

a*(*,<l>) := a(*,<f>) + ((mj x *), 0) + £C e A; (V$, V^) 

and the functional 

V&) := C e (VmiV<j>) + (H14>) + (ttK),^). 
Obviously, I^'(') * s linear, while a J (-, •) is bilinear and positive definite, since 

a\</>^) = a\ml 2 ^ + eCM\mh H - 
Hence there exists a unique £ /C 3 solving (10a). For the second equation (10b), we have 

h 

to consider the bilinear form b(-, •) : (X h , y h ) x (X h , 34) — > R defined by 

6((*,e), (V,0) := f (*,V0 - (*, v x o + <t(xuP,1>) + y(©>0 + (v x e,v) 

which is continuous and positive definite, since 

K(0,C)(0,O) = j(*>4>) - (C, v x 0) + a( x ^,0) + ^(C,C) + (v x d>,{) 
= |(^,0) + <7( X ^,0) + f (CO 



and the functional 



^IWIi-W + f HCII^n) + "ll*llk< w) 



tf'(W>0) :=-(J J ,^)-/ioK,C) 



which is obviously linear. Due to finite dimension, there is a unique solution (Ej +1 ,Hj +1 ) 
of (10b). As in Lemma 4, we see that step (Hi) of Algorithm 3 is also well-defined. □ 



5. Main theorem & Convergence analysis 

In this section, we aim to show that the two preceeding algorithms indeed define convergent 
schemes. We first consider Algorithm 3. 

5.1. Main result. We start by collecting some general assumptions. Throughout, we 
assume that the spatial meshes Th\u> are uniformly shape regular and satisfy the angle con- 
dition 

/ VO • VCj < for all hat functions Q, Q G S 1 (T h \ u ) with i ^ j. (11) 

Juj 

For x G f2 and t G [tj,tj+i), we now define for 7^ G {m^, H£, E^, J e , v e h } the time approxi- 
mations 

7M(t,x) :=^7r i (x) + ^-^ 7 ^(x) 

+ i, \ , u ^ (12) 

7«(*,x) := 7i(x), 7, + fc (t,x) := tJ +1 (x), 7a*(*,x) := 7 £ +1/2 (x) 



7l +1 (x)+7i(x) 



2 

We suppose that the general energy contribution tt(-) is uniformly bounded in L 2 (w T ), i.e. 

\W*)\\h M < C„ (13) 
with an (h, /c)-independent constant C n > for all n G L 2 (uj t ) with HnU^^) < 1 as well as 
7r(n/jfc) — 7r(n) weakly subconvergent in L 2 (ojt) (14) 

provided that the sequence n hk n is weakly subconvergent in H 1 (w r ) towards some 
n G H 1 (w T )- For the initial data, we assume 

— ^ m° weakly in L 2 (w), (15) 

as well as 

H°^H° and E° E° weakly in L 2 (fi) (16) 
Finally, for the field J, we assume sufficient regularity, e.g. J G C([0, T]; L 2 (f2)), such that 

J ± ^ J weakly in L 2 {Q T ). (17) 

Remark. Before proceeding to the actual proof, we would like to remark on the before 
mentioned assumptions. 

(i) We emphasize that all energy contributions mentioned in the introduction fulfill the 
assumptions (13)-(14) on tt(-), cf. [BSFFGPPT2]. 

(ii) As in [BSFFGPP'12], the analysis can be extended to include approximations Hh of the 
general field contribution it. In this case, one needs to ensure uniform boundedness of 
those approximations as well as the subconvergence property ^(n^) — ^ 7i~(n) weakly 
in L 2 (ujt) provided n hk is weakly subconvergent to n in H 1 (cu T ). 

(iii) The angle condition (11) is a somewhat technical but crucial ingredient for the con- 
vergence analysis. Starting from the energy decay relation 



Vml 



it has first been shown in [B'05],£/ia£ (11) and nodewise projection ensures energy 
decay even on a discrete level, i.e. 

V/ f m .)\ 2 < I |VJ„ml 2 . 



m 



This yields the inequality || Vm^ +1 W^^) — ll^ m i + ^ v /iIIl 2 (u)> which is needed in the 
upcoming proof. 

(iv) Note that assumption (11) is automatically fulfilled for tetrahedral meshes with dihe- 
dral angles that are smaller than ir/2. If the condition is satisfied by %, it can be 
ensured for the refined meshes as well, provided e.g. the strategy from [V'96, Section 
4.1] is used for refinement. 

□ 

The next statement is the main theorem of this work. 

Theorem 6 (Convergence theorem). Let (m^, v^, H^, E^) be the quantities obtained by 
either Algorithm 2 or 3 and assume (11) -(17) and 6 e (1/2,1]. Then, as (h,k) — > (0,0) 
independently of each other, a subsequence of (m^, H^, E^) converges weakly in H 1 (wqp) x 
L 2 (Qt) x L 2 (f2^) to a weak solution (m, H, E) of MLLG. In particular, each accumulation 
point o/ (m^, H/jfc, E^) is a weak solution of MLLG in the sense of Definition 1. 

The proof will roughly be done in three steps for either algorithm: 

(i) Boundedness of the discrete quantities and energies. 

(ii) Existence of weakly convergent subsequences. 

(iii) Identification of the limits as weak solutions of MLLG. 

Throughout the proof, we will apply the following discrete version of Gronwall's inequality. 

Lemma 7 (Gronwall). Let /c , . • . , > and ao, . . . , a r _i, b, C > 0, and let those quan- 
tities fulfill a < b and ai < b + C*X]j=o^i a j f or ^ = 1> ■■■,r. Then, we have a t < 
C exp (C Y!j^ k j) for£ = l,...,r. 

5.2. Analysis of Algorithm 3. As mentioned before, we first show the desired bounded- 
ness. 

Lemma 8. There exists k > such that for all k < k , the discrete quantities (m J h , E J h , H 3 h ) e 
M h x X h x y h fulfill 

3-1 3-1 

+ KIIl» + WKWlw + WKWb^ + (e- i/2)k 2 £ ||Vv* |£ 2(u;) 

i=0 i=o (lg) 



+ £ (l|H^ +1 - H* 11^^ + p* 1 - El\\l {Q) ) < C 3 



i=0 



for each j — 0, . . . , N and some constant C 3 > that only depends on on as well as 
on C n . 
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Proof. For Maxwell's equations, i.e. step (iii) of Algorithm 3, we choose (^> ft , £ ft ) = (E^ +1 , H^ +1 ) 
as special pair of test functions and get from (10b)-(10c) 

|(E^ +1 - El E^ +1 ) - (H* 1 , V x E^ +1 ) + a( X „K + \ E^ 1 ) = -(J*, E^ +1 ) and 
^(H^ 1 -Hl,H^) + (V x E*+\H^) = -^(vi.Hf). 

Summing up those two equations (and multiplying by 1/C e ), we therefore see 



£ ° ^TP»+l T?« T? i + 1 \ _i_ ° ll"tT' i + 1 ll 2 _L /'Xji+1 tt! Tji+1\ 

ZTTl^ft J + 7^11% II L»(w) + rpT^h 111 1 > 

MO/ j TTM , MO /_ i TT i TTZ+1\ 1 



(19) 

v£,H£) + ^(vi,IT h - H^ 1 ) - 7r (J i ,Ej+ 1 ). 



C\ ft' ft/ 1 V ft? ft ft / 

e 



The LLG equation (10a) is now tested with y>j = v l h e /C m ; . We get 

«K, vi) + ((ml x v *), v*) = -C e (V(m* + 6>fcv£), Vv ft ) + (H^, v ft ) + ( w (m*), v*), 

V v ' 

=0 

whence 

^KIIL H + ^ 2 ||wl||^ H = -MVml, wi) + A ( ir fc , v *) + ^(ttKJx). 

Next, we follow the lines of [A'08, ATI, BSFFGPPT2] and use the fact that || Vm^ +1 || 2 2(w) < 
||V(m^ + fcvJJH^u) stemming from the mesh condition (11), cf. [B'05], to see 

< ^||Vml||^ M + ^(Vml,Vvl) + ^||Vvl|| L2H 
= ^llVmiH^) -(9- l/2)A; 2 ||Vv^||^ H (20) 
" ^IIVftllLn + ^(HL vl) + ^(7r(m l j,v ft ). 
Multiplying the last estimate by //oA an d adding (19), we obtain 

|(l|Vm^|| 2 2H - UVmill^) + (6- l/2)^k\\Wi\\l {ui) + ^II^IlL^ 

< ^(H ft - H*+\ v*J - i-(J*,Ej+ 1 ) + ^(ttKJ, vi). 

L/ fi L/c L/f 



Next, we recall Abel's summation by parts, i.e. for arbitrary ^ el and j > 0, there holds 

1, ,2 1 1 

2 W + 2 



J 111 . J , 

53 («i - «i-l,Ui) = ^|Mj| 2 - ^|Mo| 2 + 9 53 ' Mi _ M i-l! 2 - ( 21 ) 



=1 
11 



Multiplying the above equation by k, summing up over the time intervals, and exploiting 
Abel's summation for the E ft and scalar-products, this yields 

f IIVmiH^) + (6- 1/2)M 2 £ l|Vv^||^ H + ^ g H\\b( u) + £\\K\\h m 

j—i , j—i j—i 

+ 2C~ ^ ~ -^ftllL^n) + 7T / j W^h Hl2( w ) + ^"11^1^2(0) + ^- 2^ 11 ^ -Jrl ftllL2(n) 

e i=0 e i=0 e e i=0 

< ^ £(H, - H*\ v* ) - A £(J\ E^ +1 ) + £ (7r(m*), v* ) 



^ / j\ n ft ' ft/ /-( / ^ v" 'ft / 1 /-< 

e i=0 

£ ° 11^0 112 , ^0 

2(j ll^ftllL 2 ^) 



+ ^ l|Vm°||L2 (w ) + ^°rl|E°||L2 (n ) + ^rl|H°J|L2 (n) 



v 

for any j G 1, . . . , N. By use of the inequalities of Young and Holder, the first part of the 
right-hand side can be estimated by 

^ E(H ft - H*+\ v£) - A g(J*,E^) + ^ 5>(m<), v<) 

e i=0 e i=0 e i=0 

^(IkKjll^H + I|h; +1 - ir fc ||^ (n) ) + ^ £ 

e i=0 e i=0 

^ ^ ^ 1 1 Tj 1 ^ - ! - ! 1 1 2 Z^A/ ^ > I. -m-^ 1 1 2 

e i=0 e i=0 

for any e, z/ > 0. The combination of the last two estimates yields 



j'-i , j'-i 



/l ° llVmJll^^ + (0 - 1/2)^^ £ l|Vvi||L, (w) + ^ £ ||v ft ||L H + ^-||Ei||^ (n) 

i=0 e i=0 e 

_L £ ° X^IIP 1 ^ 1 F li II 2 _L \^HT7 i + 1 ll 2 I ^0 mill 2 i ^° \^||tp+l tltj II 2 

+ 2C~ 2^ II ft - ^ftllL 2 ^) + 2^ II ft IIl2( w ) + g^ - " "ftHL 2 ^) + 2C r ^ h _ ft H L2 ( n ) 



^ ^ E(ll-Kft)||L H + l|H, +1 - H*||k (n) ) + ^ £ K 



^r^^^uin" 1 ft;iiL2( t u) -i- n iA ft - Ai ftiiL2(f7)J ' -Q- ii v ftiiL2( w ) 

e i=0 e i=0 

+ 4^5Z ll E ft +1 HL2(f7) + H Ji |lL2(n) +^°- 

e 1=0 6 i=0 

Unfortunately, the term X^~o ll^ft +1 |lL 2 (n) on ^ ne right-hand side cannot be absorbed 
by the term X^=o ll-^ft +1 IIl 2 (w) 011 the hand-side, since the latter consists only of con- 
tributions on the smaller domain u. The remedy is to artificially enlarge the first term 
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by 



7 J - 1 7 j~ 1 7 J - 1 

V^in?»+l||2 ^ fc M-pi+l T?i ||2 i fc 1 1 -pi ||2 



and absorb the first sum into the corresponding quantity on the left-hand side. With 
C==f <«-*), ^^(l-A), and * = ^(«- £), 



this yields 

% == ^llVmiH^) + (0 - l/2)M 2 i: ||Vvl||£ 2H + C v £ R||L H + ^||Ei|| 2 2(Q) 



i=0 i=0 e 

j_r< \^ Wpi+l T?i ||2 _i_ Vin?»+l||2 i ^0 IIttJ || 2 , ^ llxji+l XJM|2 

+ o E / , ll-ftfr - -^ftllL 2 ^) + 7T / j W^n IIl 2 h + 27t~" ^''l 2 ^) H 2^H ^ _Jrl /illL 2 (n) 

i=0 e i=0 e ?=o 

J-l 7 7 J-l 



e „-_q e ■ - r 



-v — 

=:b 



k J_1 

<b+ — Vtii. 

u i=0 

In order to show the desired result, we have to ensure that there are choices of e and v, 
such that the constants C v , Ch, and Ce are positive, i.e. 

(a - e) > 0, (1 - £-) > 0, and (e - -) > 

which is equivalent to ko/2 < e < a and v > ko/so- The application of the discrete of 
Gronwall inequality, from Lemma 7 yields a,j < M and thus proves the desired result. □ 

We can now conclude the existence of weakly convergent subsequences. 
Lemma 9. There exist functions (m, H, E) e H^wt,^ 2 ) x L 2 (Q t ) x L 2 (fi r ) such that 

m hk — ^ m in H 1 (a;r), (22a) 

m hk ,m± k ,m hk -± m in L 2 (H 1 (w)), (22b) 

m hk , m^,, m hk ->■ m in L 2 (w T ), (22c) 

H hk ,n± k ,H hk - Hm L 2 (fi T ), (22d) 

E feA; , E± , E/jfc - H in L 2 (fi T ), (22e) 

where the subsequences are succesively constructed, i.e. for arbitrary mesh-sizes h — > and 
timestep- sizes k — > i/iere exzs£ subindices he, kj> for which the above convergence properties 
are satisfied simultaniously . In addition, there exist some v e L 2 (u;t) with 

v hk ""vm L 2 (w T ) (23) 

for the same subsequence as above. 
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Proof. From Lemma 8, we immediately get uniform boundedness of all of those sequences. A 
compactness argument thus allows us to succesively extract weakly convergent subsequences. 
It only remains to show that the corresponding limits coincide, i.e. 

\im>y hk = lim7^ = lim7+ fe = lim^, where >y hk G {m hk , H hk , E hk }. 

In particular, Lemma 8 provides the uniform bound 

E IK 1 - m Hl^H < C3. 

Here, we used the fact that ||m^ +1 — nr^l 2 ^,-, < ^ 2 |I v 1|Il 2 M' see e -S- [A'08] or [G'12, Lemma 



3.3.2]. We rewrite 7^ G {m hk ,E hk ,U hk } as 7^ + ^{j^ 1 - 7/0 on [*j_i,fj] and thus get 

hftfc - t^IIl^^) = E / K + r J ~ ^h +1 - i h ) - 

<*EiW +1 -*(o)-^° 

and analogously 



\lhk - ltk\\l*{n T ) = E / 11^ + ~l^^h +1 - 7*) - T^lli^n) 

JV-l rt j+1 

<E / 2117^-7^11^ 



JV-l 

< 2^117^-7^11^)^0, 
3=0 

i.e. we have lim7^ fc = \im^f hk G L 2 (f2 T ) resp. L 2 (cu r ). In particular it holds that lim7 ftfc = 
lim7/ l / c . From the uniqueness of weak limits and the continuous inclusions H 1 (w-r) C 
L 2 (H 1 (w)) C L 2 (cj^), we then even conclude the convergence properties of m^,m^ fe , and 
m^fc in L 2 (H 1 (w)) as well as — ^ m in H 1 (wr). From 

HH - l|| L 2(u, T ) < |||m| - |m^||| L2(WT) + |||m^| - 1|| L > T ) 

and 

lll m wfe(*,OI - 1||l»( w ) < /imax||Vm J J| L 2 H , 

we finally deduce |m| = 1 a.e. in ojt- □ 

Lemma 10. The limit function v G L 2 (wy) equals the time derivative of m, i.e. v = <9 t m 
almost everywhere in ojt 

Proof. The proof follows the lines of [A'08] and we therefore only sketch it. The elaborated 
arguments can be found in [G'12, Lemma 3.3.12]. Using the inequality 

\\d t ra hk - v^|| L i( WT) < 2^ll v ftjL2(u; T )> 
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we exploit weak semicontinuity of the norm to see 

||<9 t m - v|| L i (WT ) < liminf ||9 t m ftfe - v^|| L i (£JT ) = as (h, k) — > (0,0), 
whence v = <9 t m almost everywhere in ojt- □ 

Proof of Theorem 6. For the LLG part of (2), we follow the lines of [A'08]. Let ip G C°°(oj t ) 
and (V>,C) e C c oo ([0,T);C oo (n)nH (curl,fi)), be arbitrary We now define test functions by 
(<f> h ,tphXh)(t, •) := (T h (m~ k x (p).l Xh ip.ly h C)( f - ')■ Recall that the L 2 -orthogonal projection 
Xy h : L 2 (Q) — > y h satisfies (u — ly h u,y h ) = for all G 34 and all u G L 2 (Q). With the 
notation (12), Equation (10a) of Algorithm 3 implies 

«/ [ (Kxvj^^-Cef (V(m^ + 0fcv^),V^)) 

Jo Jo Jo 

+ / / (7r(m^),^) 

With <f> h (t,-) := Xh( m hk x V 3 ) (A ') an( i the approximation properties of the nodal interpolation 
operator, this yields 

/ ((ovw + »i/,A- x v^), (m^ x<p))+k9 [ (Vv U) V(m^ x V )) 
Jo Jo 

+ C e / (Vra u , V(m M xtp))- I (H M „ (m ftfc x <p)) - [ (Tr(m^), (m^ x p)) 
Jo Jo Jo 

= 0(/i) 

Passing to the limit and using the strong L 2 (wT)-convergence of (m^ fc xy>) towards (mxy>), 
we get 

/ ((" v ^ + m wfc x ^hk)A m hk X¥»)) — > / (i am t + m x m t ),(m x <p)), 
Jo Jo 

kO [ (Vv^, V(m- x y>)) — ► 0, and 
Jo 

/ (Vm Ul V(m u x^))^f (Vm,V(mxv>)), 
Jo Jo 

cf. [A'08]. For the second limit, we have used the boundedness of ^llVv^^^ for 6 G 

(1/2, 1], see Lemma 8. The weak convergence properties of H^ fc and 7r(m^.) from (14) now 
yield 

/ ( H /,A- (™-hk x ¥>)) — > / (H, (m x v?)) and 
Jo Jo 

/ W%),Kxy)) — >■ / (7r(m),(mx v?)). 
Jo Jo 



'0 
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So far, we thus have proved 

/ ((am ( + m x m t ), (m x (p)) = —C e / (Vm, V(m x cp)) 
Jo Jo 

+ [ (H, (m x <p)) + I (7r(m),(m x <p)) 
Jo Jo 

Finally, we use the technical results 

(m x in, ) • (m x ip) = m t • ip, 

m, • (m x <p) = — (m x m t ) ■ ip, and 

Vm x V(m x <p) = Vm • (m x V(f) 

for the left-hand side, resp. the first term on the right-hand side to conclude (2). The equality 
m(0, •) = m° in the trace sence follows from the weak convergence m^ — ^ m in H 1 (w-r) and 
thus weak convergence of the traces. Using the weak convergence mj — ^ m° in L 2 (u), we 
finally identify the sought limit. For the Maxwell part (3)-(4) of Definition 1, we proceed as 
in [BBP'08]. Given the above definition of the testfunctions, (10b) implies 

rj~i rj~\ rj~\ rj~\ 

so [ ((E hk ) t ,^ h ) - [ (H+,Vx^)+<t/ (* w E+,iM = [ {J^+h) 
Jo Jo Jo Jo 

to I ((U hk ) t Xh)+ [ (VxE+,a) = -^ / (v^,Cfc). 
Jo Jo Jo 

We now consider each of those two terms separately. For the first term of the first equation, 

we integrate by parts in time and get 

f ({E hk ) t ^ h ) = - f (E hk , (4> h ) t ) + (E hk (T, -),MT, •)) ~(KM^ •)) 
Jo Jo v v ' 

=0 

Passing to the limit on the right-hand side, we see 

[ T ((E hk ) t ,4> h ) — ► - / T (E,^)-(E°,^(0,-)), (24) 
Jo Jo 

where we have used the assumed convergence of the initial data. For the first term in the 

second equation we proceed analogously. The convergence of the terms 



f (Hi.Vx^)-) Ah,Vx^), 
Jo Jo 

[ (x w E+,^)— ► / (X.E,V), 
Jo Jo 

[ (Jhk,4>h) — »• / ( J >^)> and 

Jo Jo 
Jo Jo 



is straightforward. Here, we have used the approximation properties (6)-(7) of the interpo- 
lation operators for the last two limits. It remains to analyze the second term in the second 
equation. Using V x E^ k (t) e 34 and the orthogonality properties of Iy h , we deduce 
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/V x E+ Xh) = /V x E+ ,C) - f (V x E+ , (1 -X y JC) 
jo Jo Jo 

= / T (VxE+,C)= Ae+,VxC)^ Ae,VxC). 
Jo Jo Jo 

For the last equality, we have used the boundary condition ( x n = on <9f2T and inte- 
gration by parts. This yields (3) and (4). 

It remains to show the energy estimate (5). From the discrete energy estimate (18), we 
get for any G [0,T] with f G [tj,t j+1 ) 

l|Vm+ (0||L H + KJL K ,) + ||H+ (f)||^ (n) + WKkWWhw 

= l|Vm+ (OII^h + f l|v,- fe ( S )||^ H + ||H+ (011^(0) + \\K k {t')\\l^) 

Jo 

< l|Vm+ {t')\\l^ + T +1 ||v, fe ( S )||^ 2H + ||H+ {t')\\l^ + ||E+ (0||L. (n) 

Jo 

Integration in time thus yields for any measurable set 3 C [0,T] 

jf ||Vm+ (f)||k (w) + jf KJL K ,) + jf ||H+ (011^^) + jf \\K k (t')\\h m <jc, 

whence weak lower semi-continuity leads to 

J II^HIlV) + J \\ m t\\l2 {uJtl) + J ||H||^2 (n) J ||E||£ a(n) < Jc 3 . 

The desired result now follows from standard measure theory, see e.g. [E'09, IV, Thm. 
4.4]. ' □ 



5.3. Analysis of Algorithm 2. This section deals with Algorithm 2 and the analysis 
follows the lines of Section 5.2. As before, we first need boundedness of the involved discrete 
quantities. 

Lemma 11. The discrete quantities (m^,E^,H^) G Mh x X h x y h fulfill 

3-1 3-1 

\\V<\\b( u) + kJ2 HWW) + \\K\\l(n) + \\K\\h(n } + V2)A: 2 £ ||V<|£ 2(u;) < C 4 

i=0 i=0 

(25) 

for each j — 0, . . . , N and some constant C4 > that depends only on \ fl\, and C n . 

Note, that in contrast to Lemma 8 from the analysis of Algorithm 3, we do not have 
boundedness of E^dW 1 - H*J|* 2(n) + ||E^ +1 - E'J^n)) in this case. 
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Proof. As before, the proof relies on the choice of the correct test functions. For Maxwell's 
equations (9b)-(9c), we choose (tphXh) = (E^ +1 ^ 2 , H^ +1//2 ) and obtain after summing up 

A ( £ ° IIT?^ 1 !! 2 _l ^° HM^+ill 2 "\ _l /tIU/ f? i+1 / 2 \\ 2 

a t{-^\\ tL 'h \\L 2 (n) + "y" 71 " L2 ( n )'' 71 Hl 2 (^) 

= -(j^^K +1/2 )-M<K +1/2 Y 

For the LLG equation (9a), we again test with <p h = V l h e /C m » and argue as in (20) to see 
gOlVm^ll^ _ ||Vml||^ H ) + ^Hvill^ + (6- l/2)lM>k\\W h \\l M 

<^(H? 1/2 ,vi) + ^(7rK),v*). 
The combination of the last two estimates thus yields for any e, v > 

j iiyy i+l II 2 , £ iipi+1112 , ^0 ||TTi+l||2 \ , G ||_, xr.i+l/2|.2 

"Hyll vm^ ||l2( w ) + 2C"" ft " L2 ( n ) 2(7~ fe IIl 2 (q)J + ^rllX^h llL 2 (n) 

+ (61 - l/2)Ml|Vv £11*^ + g(a - e)||vi||L, (w) 

^ ^ || T i+l/2||2 I 1 llP^ +1 / 2 ll 2 I ^° H^mMI! 2 

Multiplying by /c and summing over the timesteps, we see 



^"Vrr^'ll 2 _l_ £ ° llT^'ll 2 4- ^° lll^'ll 2 -L^X^IKx F i+1 / 2 H 2 

vm ^llL 2 H + ^II^Hl 2 ^) + ^-II^IIl 2 ^) + ^2j llXw-^h IIl 2 (^) 



We i=0 



+ (0 - 1/2)^^ £ ||Vv* + ^*(a - e) £ K||2, (w) 



i=0 e i=0 

<T ^ V II T i+1 / 2 || 2 I fc 1IF' +1/2 II 2 I ^0 IUrrr.MII 2 

e i=0 e i=0 e i=0 

_L ^° HY7rv,0 II 2 1 £ ° IITPOII 2 1 ^° lluOll 2 

-t- — 1| vm ft || L 2( w ) -t- ^tII^Hl 2 ^) "+" 2^7~" n /illL 2 (n) 

~ C+ 8^C~ e ^ (ll E ^ +1 HL 2 (n) + Philip))- 

6 i=0 

Analogously to the last section, the term X^=o II E'/i +1/ ' 2 II l 2 (^) * n the third line cannot be 
absorbed by the one on the left-hand side directly due to the different domains. We thus 
again extend the quantity by 

j'-i j'-i 

(H E /i +1 |lL 2 (n) + II E IIIl 2 (Q)) < ll E jllL 2 (m + ll E fellL 2 (ffl> 
i=0 i=0 

and absorb the first part by the left-hand side for appropriate v. As before, the assertion 
then follows by an application of the discrete Gronwall's inequality. □ 
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Analogously to Lemma 9, we conclude the existence of weakly convergent subsequences 
that fulfill 



m hk m in H (u T ), 


(26a) 


m hk ,m± k ,m hk — m in L^H 1 ), 


(26b) 


m hk ,m± k ,m hk -> m in L 2 (u T ), 


(26c) 


H^fc, H^ fe , H hk — ^ H in L 2 (f2 T ), 


(26d) 


Em, E^ fc , Em — ^ E in L 2 (r2r), 


(26e) 


v m v in L 2 (w r ). 


(26f) 



Note also that the above mentioned boundedness of 

j-i 

El'llTT^ 1 "H"* II 2 _l_ IIT?^ 1 T?* II 2 

i=0 

is not necessary to prove that the limits of the in time piecewise constant and piecewise 
affine approximations coincide. The remedy is a clever use of the midpoint rule. The proof 
of Theorem 6 for Algorithm 2 then completely follows the lines of the one for Algorithm 3. 

Proof of the convergence Theorem 6 for Algorithm 2. Again using (0^,^, Ch)(t, •) := (X/j(m^ fc 

<p),I Xh tl>,Iy h C)(t,-) for any (<^,C) e C°°(u> T ) x C c °° ([0, T); C°°(n) n H (curl, fi)) 2 , Algo- 
rithm 2 implies 



x 



(27) 



VmA) + / ((m ftfc x v hk ),<f> h ) = -C e I (V(m^ + ekv hk ),V<f> h ) 
Jo •/ o 

+ / (H hk ,<i> h )+ [ (n(m hk ),(j) h ) 
Jo Jo 

e I ((E hk ) t ,if> h ) - [ (Hhk, V x ij) h ) +<j / = - / (3hk,iph) 

Jo Jo Jo Jo 

Ho [ ((H hk ) t ,Ch) + / (V x E^CJ = -/i / (Vft fc ,Ch)- 
Jo Jo Jo 

Next, we use the fact that the above scalar products containing E hk and H hk can be expressed 
by means of E hk and H^, respectively, by use of piecewise constant test functions in time. 
For A G C°°(Qt), consider the piecewise constant approximation A~ G V°(I k , C°°(Q)) with 
A~(t) = A(tj) for i G [tj,tj + i). Since the midpoint rule is exact for the (piecewise) affine 
function (H^, A~), there holds 

N-l rt , +x i N-l / 



/ (H hk ,A-) = J2 ^(Hf 1 + H^,A(t J ))=A;^(H ?ifc ,A-) 

Jo j=0 A ^ j=0 V 

= J2 / (H ftfc ,A-) = / (H^A-). 

-_n Jtn JO 



3=0 " l i 
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Analogously, we get 

I (H U) Vx A-)= / (HfcVxA"), 
Jo Jo 

[ ( x jE hk ,A-)= [ (x w Ewfe,A-), 
Jo Jo 

(VxE u ,A> / (VxE^A-). 



Now, let (j) h ,i/) h Xh e ^ ,0 (^fc) denote the in time piecewise constant approximations of <f>h, if>h, 
and £h respectively. We then get 

/ (H hk ,<j> h ) = / (H hk ,fa) + / (ii hk ,(j) h - fa) 
Jo Jo Jo 

= / (H ftfc ,$~) + / (H hk ,<f> h -fa) 
Jo Jo 

" v ' 

= '- a hk 

In the following, we are going to show lim(/ ljfc )_ 5> ( 0i0 ) a hk = 0. By definition, we get 

a hk = [ (H hk ,I h (m- k x<p- m hk x <p~)) 
Jo 

= / (Hhfc, (m- x v? - m- fc x p")) + / (H fcfc , (1 - X h )(m^ x <p - x yT)). 
Jo Jo 

For the second term, we immediately get 

/ (H hk , (1 - T h ){m- hk x<p- m- x = G(/ i 2 ) 
Jo 

due to boundedness of ||H hfc || L 2( n ^, ||m^(-) || H i( n ) and ||<p|| W 2,<x>, and therefore elementwise 
boundedness of ||(m^ fe x (p) — (m^ k x ^~)||h 2 (t) fc> r an Y T <E %■ For the first term on the 
right-hand side, the mean value theorem yields 

rtj+i rtj+i 
/ \\( m hk >< V>) - ( m hk X <P~)\\h(fl) < \\<P~ ~ ¥>llL~(fe,t J+1 ];L2(fi)) / H^Hl^O) 

2 2 /^-H 

whence strong L 2 (fi T )-convergence of (m^ x <p — m hk x <p~) to 0. Altogether, we thus get 

/ (Hwfc, (m^ x v? - m fefc x p")) < ||H ftfc || L2(nT )||(m^ x <p) - (m^ x yT)|| La( n T ) — ► 0. 
Jo 

Analogously, we derive 0^ — >■ m x <p strongly in L 2 (Qt) and thus conclude 

/ (HhkM—tf (H.m x v?) as (M)^ (0,0). 
Jo Jo 
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Similar arguments show 

rj~\ rj~\ rj~i 

I (H U) Vxi)=f (H ( ,.Vx^)+ /' (Hft fc) Vx(^-^))^/ (H.Vxt/O. 
Jo Jo Jo Jo 

I {xJ^hk^h) = [ (xAk,il>h)+ I (xAk,tl>h-*l>h) — ►/ (x«E,^), 
Jo Jo Jo Jo 

[ (VxE hk ,Ck)=[ (VxE^Q+f (VxE^C-C")^ / (E,Vx(). 
Jo Jo Jo Jo 

Using the convergence properties (26a), the remainder of the proof follows as for the one of 
Theorem 6 for Algorithm 3. As for the energy estimate, we again utilize weak semi-continuity 
and the discrete energy estimate (25). □ 

6. Numerical examples 

We study the standard /i-mag benchmark problem no. 4, see [muM AG] using Algorithm 2 
and Algorithm 3. Here, the effective field consists of the magnetic field H from Maxwell's 
equation and some constant external field H ex t, i-e. 7r(m^) = H ex t for all j = 1, ... , N. This 
problem has been solved previously using the midpoint scheme in [B'10], and we also use 
those results for comparison. 

Despite the fact that the system (9) in Algorithm 2 is linear, for computational reasons 
it is preferable to solve Maxwell's and LLG equations separately. After decoupling, the 
corresponding linear systems can be solved using dedicated linear solvers. This leads to a 
considerable improvement in computational performance, cf. [BBP'08]. In order to decouple 
the respective equations in (9), we employ a simple block Gauss-Seidel algorithm. For 
simplicity we set o = 0, J = 0. Assuming the solution v^ -1 , H^, is known for a fixed 
time level j, we set G° = H{, F° = E{, and w° = v^ 1 and iterate the following problem 
over I: Find w e h , F e h , G{ e K* x X h x y h such that for all 4> h ^hXh G K.* x x h x Vh, we 

h h 

have 



a(wl<f> h ) + ((mi x wj),* fc ) = -C e (V(mJ + 9kw[), V<f> h ) 

+ (Gl 1 + ii ext ,<i> h ) 



(28a) 



eoliKM - (Gl VxW = eo\{KM (28b) 

IXo\{GUh) + (V x F{,Ch) = ^l^uXh) ~ »o{<Xk) (28c) 
t-l rit t-it-l -ci TTi-e— ll 



until |K - w^loo + \\G f h - Grioc + \\K - F^loo < TOL. In this setting, F* is an 
appr 
have 



approximation of E^ +1 ^ 2 and G e h is an approximation of H{ +1 ^ 2 , respectively. Therefore, we 



9 9 T?3 +1 _ TfJ 

2 -(Fi - E{) « ^(Ef 1 / 2 - E{) = = d t E{ + \ 

Analogous treatment of the H{ +1//2 -term thus motivates the above algorithm. We obtain the 
solution on the time level j + 1 as v j h = w{, H( +1 = 2G f h - H{, E{ +1 = 2F e h - E{. The 
linear system (28a) is solved using a direct solver, where the constraint on the space JC m j 

h 
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is realized via a Lagrange multiplier, see [GHMPS]. For the solution of the linear system, 
(28b)-(28c) we employ a multigrid preconditioned Uzawa algorithm from [BBP'08]. 

The physical parameters for the computation were ji Q = 1.25667 x 1CT 6 , e = 0.88422 x 
10" 11 , A = 1.3X10- 11 , M s = 8xl0 5 ,7 = 2.211xl0 5 , a = 0.02, H ext = (/i M s )- 1 (-24.6, 4.3, 0), 
C e = 2A(/i M^)~ 1 . Here, 7 denotes the gyromagnetic ratio, and M s is the so-called satu- 
ration magnetization, see e.g. [BSFFGPP'12]. We set 9 = 1 in both, Algorithm 2 and 3. 
The ferromagnetic domain u = 0.5 x 0.125 x 0.003 (/zm) is uniformly partitioned into cubes 
with dimensions of (3.90625 x 3.90625 x 3) (ran), each cube consisting of six tetrahedra. The 
Maxwell's equations are solved on the domain Q — (4 x 4 x 3.072) (fim). The finite ele- 
ment mesh for the domain Q is constructed by gradual refinement towards the ferromagnetic 
domain u, see Figure 1. We take a uniform timestep k = 0.05 which is two times larger 
than the time-step required for the midpoint scheme [B'10]. Note that the scheme admits 
time-steps up to k — 1, the smaller time-step has been chosen to attain the desired accuracy. 




Figure 1. Mesh for the domain Q at £3 = (left) and zoom at the mesh for 
the domain u> at x 3 = (right). 

The initial condition m for the magnetization is an equilibrium "S-state", see Figure 2, 
which is computed from a long-time simulation as in [BBP'08], [B'10]. The initial condition 
H is obtained from the magnetostatic approximation of Maxwell's equation with and E = 
0, for details see [B'10]. In Figure 3 we plot the evolution of the average components m x and 
ni2 of the magnetization for Algorithm 2 and Algorithm 3. For comparison, we also present 
the results computed with the midpoint scheme from [B'10] with timestep k = 0.02. 

We also show a snapshot of the magnetization for Algorithm 2 and the midpoint scheme 
at times when J^m\{t) = in Figures 4 and 5, respectively. We conclude that the 
results for both algorithms are in good agreement with those computed with the midpoint 
scheme. 
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Figure 2. Initial condition m°. 
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Figure 3. Evolution of \u\ 1 f mi and \u\ 1 J w m2, where rrij denotes the 
j-th component of the computed magnetization m : u — > K 3 . 
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